DOCENTES

Contenido del curso

• T1. Regresión Lineal, suposiciones y requisitos.
• T2. Estimación de los parámetros e interpretación de los valores.
• T3. Pruebas de hipótesis relacionadas con los parámetros.
• T4. Evaluación de modelos de acuerdo a las suposiciones de normalidad e independencia.
• T5. Transformaciones de las variables, y sus implicaciones en las pruebas de hipótesis.
• T6. Series de Tiempo. Estacionariedad y cómo obtener una serie estacionaria a partir de una que no lo es. Correlogramas.
• T7. Modelos de Box y Cox. Estimación de parámetros.
• T8. Criterios de evaluación de un modelo de series de tiempo. Transformaciones.

Bibliografía.

Regression Modeling and Data Analysis with Applications in R. Samprit Chatterjee, Jeffrey S. Simonoff
Montgomery, D. C., Peck, E. A., & Vining, G. G. (2002). Introducción al análisis de regresión lineal (3ª ed.). Limusa Wiley.
Time Series Analysis and Its Applications: With R Examples by Shumway and Stoffer.
Time Series Analysis: Forecasting and Control by Author(s): George E. P. Box, Gwilym M. Jenkins, Gregory C. Reinsel, Greta M. Ljung.
https://fhernanb.github.io/libro_regresion/rls.html#modelo-estad%C3%ADstico
https://bookdown.org/victor_morales/SeriesdeTiempo/

Motivación

Los métodos de regresión son técnicas estadísticas utilizadas para modelar y comprender la relación entre variables. La motivación detrás de estos métodos radica en la necesidad de analizar y cuantificar cómo una o más variables predictoras influyen en una variable respuesta. Los métodos de regresión son ampliamente utilizados en diversos campos, incluyendo la investigación científica, la economía, la medicina, la ingeniería y más. A continuación, se detallan algunas de las principales motivaciones detrás de estos métodos:

Modelado de Relaciones:

La motivación fundamental de la regresión es modelar la relación entre variables. En muchos casos, existe la necesidad de entender cómo una variable dependiente cambia en función de una o más variables independientes. Por ejemplo, en la economía, podría ser necesario entender cómo las tasas de interés afectan el gasto del consumidor.

Predicción y Estimación: Los modelos de regresión también se utilizan para hacer predicciones y estimaciones. Dado un conjunto de datos históricos, los modelos de regresión pueden ayudar a predecir valores futuros de la variable dependiente en función de los valores de las variables independientes. Esto es especialmente útil en campos como la meteorología, la finanzas y el marketing.

Control y Optimización: En algunos casos, se utilizan modelos de regresión para optimizar y controlar procesos. Por ejemplo, en la manufactura, los modelos de regresión pueden usarse para identificar las condiciones ideales que conducen a la máxima eficiencia o calidad del producto.

Validación de Teorías: En la investigación científica y social, los modelos de regresión pueden usarse para validar o refutar teorías existentes. Al comparar los resultados del modelo con las expectativas teóricas, se puede evaluar la validez de las hipótesis.

Gestión de Riesgos: Los modelos de regresión también se utilizan para evaluar riesgos y tomar decisiones informadas. En finanzas, por ejemplo, se pueden utilizar para evaluar el riesgo de inversión y la exposición a diferentes factores del mercado.

Ejemplo

El modelo CAPM (Capital Asset Pricing Model) es un modelo de valoración de activos financieros desarrollado por William Sharpe que permite estimar su rentabilidad esperada en función del riesgo sistemático. Su desarrollo está basado en diversas formulaciones de Harry Markowitz sobre la diversificación y la teoría moderna de portafolios. Dos puntos claves son:

1. El modelo CAPM es utilizado para calcular la rentabilidad que un inversionista debe exigir al realizar una inversión en un activo financiero en función del riesgo que está asumiendo.
2. El modelo CAPM establece una relación lineal entre el rendimiento esperado de un activo y su riesgo sistemático, medido por su \(\beta_i\).

La fórmula del modelo CAPM es la siguiente:

\(E(r_i)=r_f+\beta_i*(E(r_m)-r_f)\)

Donde:

\(E(r_i)\) es la tasa de rentabilidad esperada de un activo concreto. \(rf\) es la rentabilidad del activo sin riesgo. \(\beta_i\) es la medida de la sensibilidad del activo respecto a su benchmark. \(E(r_m)\) es la tasa de rentabilidad esperada del mercado en que cotiza el activo.

Análisis de Causa y Efecto: Los modelos de regresión pueden ayudar a establecer relaciones causa-efecto entre variables. Esto es útil para comprender cómo los cambios en una variable influyen en otras variables y viceversa.

En resumen, la motivación detrás de los métodos de regresión radica en la necesidad de entender, modelar, predecir y cuantificar las relaciones entre variables en una variedad de campos. Estos métodos permiten tomar decisiones informadas, realizar análisis profundos y obtener información valiosa a partir de los datos disponibles.

1.1 Lectura de datos

A continuación, cargaremos dos tablas de datos que contienen información sobre el precio del dólar, también conocida como tasa representativa del mercado, y otra con la información del precio del aguacate Hass tomada de una fuente gubernamental y de la página del SIPSA.

Los enlaces para acceder a información relacionada son:

https://www.agronet.gov.co/estadistica/Paginas/home.aspx?cod=11 https://www.bde.es/webbe/es/estadisticas/temas/tipos-cambio.html

# Cargue y limpieza de datos.

dolar <- read.csv2("datasets/Tasa_de_Cambio_Representativa_del__Mercado_-Historico.csv") %>% 
  clean_names()
hass <- read.xlsx("datasets/Hass_Precios_Historicos.xlsx") %>% clean_names()

1.2 Creación de funciones.

Se programa una función de flextable que mejorará la impresión de tablas de resumen. Para usarla sólo bastará invocar la función ftable() al final de la sentencia (usando %>% tidyverse) o con el ftable() encapsular o encerrar la tabla que se quiere ajustar.

# Funciones

## Función para crear flextable
ftable <- function(x) {
  x %>% 
    flextable() %>% 
    theme_vanilla() %>%
    color(part = "footer", color = "#666666") %>%
    color( part = "header", color = "#FFFFFF") %>%
    bg( part = "header", bg = "#2c7fb8") %>%
    fontsize(size = 11) %>%
    font(fontname = 'Calibri') %>%
    # Ajustes de ancho y tipo de alineación de las columnas
    set_table_properties(layout = "autofit") %>% 
    # width(j=1, width = 3) %>%
    align(i = NULL, j = c(2:ncol(x)), align = "right", part = "all")
}

Este es un ejemplo del uso de la función:

Sin flextable

dolar %>% 
  head(5) 

##      valor unidad vigenciadesde vigenciahasta
## 1 4,761.64    COP      22/12/22      22/12/22
## 2 4,781.28    COP      20/12/22      20/12/22
## 3 4,802.48    COP      17/12/22      19/12/22
## 4 4,836.24    COP      13/12/22      13/12/22
## 5 4,815.99    COP      10/12/22      12/12/22

Con flextable

dolar %>% 
  head(5) %>% 
  ftable()

1.3 Depuración de los datos

# Vamos a sacar un valor promedio de cada variable por mes
#Para el dolar vamos a tomar el campo vigenciadesde

dolar <- dolar %>% 
  mutate(fecha = as.Date(vigenciadesde, format = "%d/%m/%y")) %>% 
  mutate(mes = month(fecha), anio = year(fecha)) %>% 
  mutate(valor = as.double(str_replace(valor,",",""))) %>% 
  group_by(anio,mes) %>%
  summarise(precio_dolar = mean(valor), .groups = "drop") 

dolar %>%
  datatable()

hass %>% 
  head(5) %>% 
  ftable()

# 🖇️ Se genera el vector de meses para usarlo más abajo con la función match.
meses <- c("Enero", "Febrero", "Marzo", "Abril", "Mayo",
           "Junio", "Julio", "Agosto", "Septiembre", "Octubre",
           "Noviembre", "Diciembre")

hass <- hass %>% 
  # Otra forma de sacar el anio.
  #mutate(anio = substring(fecha,nchar(fecha)-4,nchar(fecha))) %>% 
  #mutate(espacio = grepl(", ",fecha))
  mutate(mes = sapply(strsplit(fecha, " "), "[", 2),
         anio = sapply(strsplit(fecha, " "), "[", 5)) %>% 
  mutate(anio = as.double(anio)) %>% 
  # 🖇 Se utiliza la función match con el vector meses.
  mutate(mes = match(mes,meses)) %>% 
  group_by(anio,mes) %>%
  summarise(precio_aguacate_kg = mean(precio_kg), .groups = "drop") 

hass %>% datatable()

hass_dolar <- dolar %>% 
  right_join(hass, by = c("mes","anio")) 

hass_dolar %>% 
  head(5) %>% 
  ftable()

1.4 Diagrama de dispersión.

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
  geom_point()+ theme_light()

Como se puede notar, hay alguna relación entre ambas variables puesto que si el precio del dólar aumenta, también lo hace el precio de aguacate. La relación no es perfecta pues no vemos una curva suave que sea capaz de pasar por todos los puntos, más bien debemos reconocer que hay ciertas variaciones en el proceso. Estas variaciones pueden deberse a factores no contemplados por el modelo, después de todo no podemos esperar que el precio del aguacate dependa exclusivamente del precio del dólar.

1.5 Generando un modelo súper simplificado (modelo ingenuo)

Ahora vamos a generar el modelo más sencillo que podamos imaginar para predecir el precio del aguacate, al cual llamaremos un modelo ingenuo. La razón del apelativo es que el modelo “ingenuamente” supondrá que es posible predecir el precio que tendrá el aguacate con el promedio histórico.

La siguiente gráfica muestra la linea horizontal que simboliza el promedio histórico del precio del aguacate. También debe notar que predecir basado en esta línea es despreciar cualquier aporte de la variable precio del dólar en la predicción, es decir, es un modelo que no explota la relación (sea ésta débil o sea ésta fuerte) entre el precio del aguacate y su variable explicativa precio del dólar.

ggplot(data = hass_dolar,
       aes(x= precio_dolar, y= precio_aguacate_kg )) +
geom_point()+ 
geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg),
           color = "red")+
geom_text(data=hass_dolar,
          aes(x= 5700, y = mean(hass_dolar$precio_aguacate_kg),
              label = paste("Precio promedio histórico: $",
                            round(mean(hass_dolar$precio_aguacate_kg), 2))),
          hjust = 1.2, vjust = -0.2, color = "red"
) +
theme_light()

Observe que el modelo ingenuo no sería tan inadecuado en dos circunstancias:

• Si los datos se encuentran muy cercanos al precio histórico de manera consistente sin importar en que valor del precio del dólar me posicione. Lo que vendría a indicar un escenario en el cual moverme a lo largo de diferentes valores del precio del dólar no genera ningún patrón de distanciamiento frente al valor histórico de referencia.

• Si los datos del precio del aguacate se alejan de la línea histórica de referencia pero el precio del dólar tampoco puede seguirlos, es decir, si la relación entre precio del dólar y precio del aguacate se parece a una nube de puntos sin ningún patrón o tendencia evidente. En este caso cualquier cosa que supongamos de la relación entre precio de dólar y precio del aguacate será producto de nuestra imaginación 🤦🤦

En estos dos escenarios es mejor retener el modelo ingenuo a falta de una variable que de verdad explique lo que pasa con el precio del aguacate.

Ahora bien, ninguno de los dos escenarios anteriores parece ser el nuestro, más bien acá se nota que el precio del dólar si está acompasado con el precio del aguacate. Sin embargo antes de entrar en la búsqueda de esas relaciones es importante notar dos cosas:

• El modelo ingenuo es mi modelo de referencia, esto significa que si un modelo predictivo construido para predecir el precio del aguacate es mejor, lo tendrá que ser respecto a este modelo ingenuo.

• Es posible medir el desajuste actual de mi modelo ingenuo tomando la distancia de cada punto a la recta horizontal de promedio histórico.

1.6 Calidad de ajuste del modelo ingenuo

Sea \(e_i = y_i - \bar{y}\) las diferencia entre el dato del precio de aguacate \(y_i\) y el promedio histórico \(\bar{y}\). Como tengo muchas diferencias, será necesario definir una métrica resumen, por conveniencia elijamos la suma de los cuadrados de las diferencias, debido a que si sumamos las diferencias, éstas se me compensarán, pues diferencias negativas cancelarán las positivas, y no quiero esto. Más bien me interesa que ambas sumen a mi medida de desajuste. Llamemos a esta medida la suma de cuadrados del error:

\[SCE_{\text{ingenuo}} = \sum^{n}_{i=1}e_i²=\sum^{n}_{i=1}(y_i - \bar{y})²\]

En la definición de la cantidad quisimos usar el subíndice “ingenuo” para nombrar al SCE, esto con el fin de hacer énfasis en que tal métrica se asocia al modelo, si el modelo cambia a otro tipo de modelo el SCE cambiará también: será la distancia de cada \(y_i\) a la curva definida por el otro modelo. Veremos esto más adelante. No obstante ya podemos sacar nuestras primeras conclusiones:

1. Entre más grandes sean cada una de las distancias \(e_i\) individuales, más grande será SCE.
2. Todas las distancias contribuyen con igual importancia al SCE, excepto en lo que se refiere a su magnitud no hay porque pensar que un distanciamiento de \(u\) unidades sea más importante que otro de las mismas \(u\) unidades. Esto puede ser obvio pero no lo es cuando, por ejemplo, queremos castigar más las distancias que se dan en la zona central del gráfico que las que se dan en los extremos del gráfico.
3. El SCE es sensible a la cantidad de datos, razón por la cual entre más datos (es decir entre más grande sea n), más grande será SCE. Esto dificulta seriamente la posibilidad de comparar modelos que fueron calculados sobre una cantidad distinta de puntos.

De acuerdo con lo anterior modifiquemos un poco nuestra definición de desajuste, así:

\[MSE_{\text{ingenuo}} = \color{red}{\frac{1}{n}}\sum^{n}_{i=1}(y_i - \bar{y})²\]

Es decir, nuestra medida de desajuste será el promedio de la suma de las diferencias cuadradas de cada dato \(y_i\) respecto a su media \(\bar{y}\), dicho de otra manera la varianza de la variable \(Y\) a predecir!!!!

A continuación invitamos al estudiante a calcular la varianza de la variable a predecir, es decir, el desajuste del modelo ingenuo.

# En esta sección el estudiante desarrollará los cálculos.

# 💡 Recordemos que la varianza de los residuales respecto del modelo ingenuo, es la misma varianza de la variable a predecir.

Observe el siguiente gráfico en el que hemos pintado las distancias que hay entre cada dato \(y_i\) observado y el promedio histórico (valor predicho \(\bar{y}\) por el modelo ingenuo)

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg )) +
  geom_point()+ 
  geom_hline(yintercept = mean(hass_dolar$precio_aguacate_kg))+
  geom_segment(aes(xend=precio_dolar, yend=mean(precio_aguacate_kg)),
               col='red', lty='dashed')+
  theme_light()

1.7 Construyendo un modelo lineal simple con los datos

Este es el momento de comenzar a usar la información extra con la que contamos, esto es, el precio del dólar, el cual se convertirá en una variable predictora del precio del aguacate. En estadística no hay una única forma de hacer esto, pero lo usual es partir de especificaciones muy sencillas. Veamos:

Sea \(Y\) el precio del aguacate y sea \(X\) el precio del dólar. Investiguemos si este modelo sencillo nos funciona:

\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, {\sigma_e}^2)\]

1.8 Algunas observaciones sobre las fortalezas y limitaciones del modelo

Expliquemos qué implica la ecuación anterior haciendo ciertas anotaciones de gran relevancia práctica para entender cómo los estadísticos piensan al momento de plantear modelos:

Observación genial 1:

Para un x fijo, es decir un x dado, es posible calcular el valor de \(y\) con el modelo anterior, donde tácitamente estamos diciendo que el valor de \(y\) dependerá de \(x\) en forma exacta, excepto por una perturbación aleatoria \(e\) que representará la parte del valor de \(y\) que no puede ser capturada por el término \(\beta_{0} + \beta_{1} * x\).

¿Cuándo sucederá que el valor de \(y\) pueda ser exactamente capturado por el término \(\beta_{0} + \beta_{1} * x\)? Cuando \(y\) dependa en forma lineal exacta del valor de \(x\). Esta es una condición demasiado fuerte que rara vez ocurre en la práctica, por eso agregamos la perturbación \(e\) como una forma de modelar la incertidumbre en la determinación de \(y\) dado un valor de \(x\). Esto no soluciona todos los problemas pero ayuda a obtener un modelo relativamente más realista.

¿Cuáles son los riesgos que implica asumir a priori un comportamiento específico para el error?

Observación genial 2:

La perturbación \(e\) puede recibir varios nombres en estadística, la encontrarás nombrada como: perturbación, choque, error, o residual. Lo importante no es el nombre que reciba, sino las condiciones que vamos a suponer sobre sus valores.

En particular no vamos a exigir que \(e\) tenga valores predecibles, muy al contrario vamos a exigir que sea un valor aleatorio, precisamente porque está modelando incertidumbre. Pero esto no implica que no podamos poner ciertas condiciones sobre el tipo de aleatoriedad deseada para \(e\). En este caso supondremos que \(e\) se distribuye como una variable aleatoria normal con media \(\mu=0\) y desviación estándar \(\sigma_e\). Conviene aclarar que esto es una suposición, no implica que efectivamente los errores vayan a cumplir tal supuesto.

Observación genial 3:

¿Por qué suponemos normalidad para \(e\)? Porque si \(e\) es en efecto una especie de componente incierto en el valor de \(y\), entonces muy seguramente será la suma de muchos efectos independientes que lo están provocando. En nuestro contexto estos efectos inciertos pueden ser:

• \(F_1\): Cambios en la productividad de los cultivos.
• \(F_2\): Condiciones climatológicas.
• \(F_3\): Precios de los insumos agrícolas.
• \(F_4\): Capacidad de negociación de los productores.
• \(F_5\): Presencia o no de subsidios estatales.
• \(F_6\): Precios de productos complementarios o sustitutos.

…

• \(F_k\): precio de la gasolina.

Y un gran etcétera.

Es decir, existen innumerables factores, distintos al precio del dólar que impiden que el precio del aguacate pueda ser determinado de forma exacta sólo por observar el valor del dólar. En estadística se ha estudiado que cuando no estamos midiendo los demás factores que afectan el valor de una variable, y dejamos que estos contribuyan de forma desacoplada e independiente a dicho valor, el efecto general \(e\) que resume el efecto global de todos los factores a la vez, suele comportarse bajo la distribución normal sencillamente porque hay un teorema en matemáticas, conocido como el Teorema del Límite Central que básicamente así lo garantiza, al postular que: la suma de \(k\) efectos aleatorios independientes tiende a adoptar el comportamiento normal conforme la cantidad k de efectos crece.

Es importante destacar que el Teorema del Límite Central tiene algunas condiciones y suposiciones, como la independencia de las variables aleatorias y que la suma debe efectuarse sobre una cantidad \(k\) de ellas suficientemente grande. Éste es uno de los conceptos fundamentales en estadística y es ampliamente utilizado en la teoría y la práctica de esta disciplina.

Observación para nada genial 4 😥😢:

Excepto por el componente aleatorio \(e\), el modelo restringe la relación de \(Y\) y \(X\) al universo de relaciones lineales. Existirá una posible relación por cada par de valores \(\beta_0\) y \(\beta_1\) que decidamos elegir. Pero por más que nos esforcemos en modificarlos siempre conducirán a relaciones representadas por líneas rectas, de ahí que el modelo asuma el nombre de regresión lineal.

Si queremos capturar otras posibles dependencias de \(Y\) respecto a \(X\) podemos generalizar la relación usando \(y = f(x) + e\) donde \(f(x)\) puede ser una nueva especificación funcional con otra estructura deseada, por ejemplo un polinomio o cualquier otra función que querramos definir. Lo importante es que detrás de la definición de \(f(x)\) haya alguna justificación producto de haber analizado los datos en búsqueda de relaciones que capturen bien la dependencia.

Observación genial final:

Son los datos observados para cada valor de \(y\) y \(x\) los que nos deben informar sobre el par de valores \(\beta_0\) y \(\beta_1\) que crean la relación lineal que mejor representa nuestra nube de puntos, pero de nada sirven los valores observados si no definimos una métrica que nos informe sobre la calidad del ajuste.

1.9 Introducción a las medidas de ajuste para modelos no ingenuos

De la observación final del apartado anterior nos queda una inquietud: ¿Cómo decidimos el par de valores a asignar a los parámetros del modelo lineal si no tenemos un criterio para poder saber cuál es el mejor modelo?

Para salir de este embrollo, reconozcamos que nunca en un problema real sujeto a incertidumbre una recta pasará exactamente por todos los puntos, razón por la cual aparecerán componentes de error \(e_i\) por cada \(y_i\) y \(x_i\) observado. Vimos que esto ocurrió con el modelo ingenuo y por supuesto ocurre también para un modelo no ingenuo. La idea es entonces reducir la magnitud de estos errores al mínimo posible y el par de valores \(\beta_1\) y \(\beta_2\) que así lo logren será nuestra elección óptima de cara al objetivo de minimización del error.

Otro criterio de ajuste que se podría usar es el de maximizar la probabilidad de ocurrencia de los valores observados bajo el supuesto de que el modelo usa valores \(\beta_0\) y \(\beta_1\) previamente especificados. De esta manera si un par de valores hace menos probable haber obtenido nuestros datos observados deberá descartarse.

¿Pero cómo podemos medir la probabilidad de ocurrencia de nuestros datos observados una vez damos valores específicos para los betas? Esta es una cuestión relacionada con los métodos de estimación de parámetros, nombre con el que se conoce al procedimiento estadístico encaminado a elegir los mejores betas para un conjunto de datos observados. Si se trata de maximizar la probabilidad de ocurrencia deberíamos usar el método de estimación por máxima verosimilitud. Se recomienda al estudiante dar una mirada a la bibliografía sugerida para entender cómo opera este método. Sin embargo podemos anticipar que ambos métodos conducen, en el caso de la regresión lineal simple, a iguales estimaciones de los paramétros.

Ahora vamos a proceder de manera más intuitiva, dejando de lado el formalismo matemático, y vamos a ejemplificar un posible modelo ajustado usando valores \(\beta_0=700\) y \(\beta_1 = 1.2\). Hemos elegido estos valores por simple inspección visual así que no esperamos que sean óptimos en ningún sentido, veamos:

# Calculamos las predicciones de "y" usando el modelo
b0 = 700
b1 = 1.2
x= hass_dolar$precio_dolar
y_pred = b0 + b1*x

hass_dolar %>% 
  ggplot(aes(x= precio_dolar, y= precio_aguacate_kg ))+
  geom_point() +
  geom_line( aes( x= x, y = y_pred), col = "red")

Aunque no estamos seguros de la optimalidad de los betas elegidos, el ajuste basado en este modelo es a simple vista decente, pero no deberíamos depender de esto. Es necesario usar las herramientas que R nos ofrece para estimar los betas. Antes de pasar a su uso ciego, revisemos en qué consideraciones estadísticas se basan estos métodos. Para ello definamos primero los errores cometidos por el modelo, y más aún, la medida de desempeño que usaremos para decidir cuál es la mejor elección de los betas. Usaremos entonces la métrica MSE(Mean Squared Errors) dada por:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}e_i²\]

Por otro lado podemos decir que \(e_i = y_i - \hat{y}\), donde la cantidad \(\hat{y}\) representa el valor de la variable \(Y\) que arroja el modelo, el cual es por regla general diferente al verdadero valor observado de \(Y\), es decir, los errores son las diferencia de lo observado \(y_i\) y lo predicho por el modelo basado en la recta \(\hat{y_i}\), por lo tanto:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \hat{y_i})²\]

Y más específicamente:

\[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - (\beta_0 + \beta_1*x_i))²\] \[MSE_{\text{NO ingenuo}} = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\]

Con esta expresión ya se hace patente que \(MSE_{\text{NO ingenuo}}\) es función de los pares de valores \((x_i, y_i)\) a los que se debe ajustar la recta, así como de los valores que elijamos para los betas.

Sin embargo, como los pares \((x_i, y_i)\) ya vienen dados por nuestra base de datos o, hablando en términos más generales, por la información recaudada para la construcción del modelo, sólo tendremos opción de variar los betas en aras de minimizar la cantidad \(MSE_{\text{NO ingenuo}}\). Dicho de otra manera los pares de valores \((x_i, y_i)\) se comportan como constantes a lo largo de mi proceso de minimización.

Dado lo anterior, es buena idea definir una función MSE en R que permita calcular, para diferentes elecciones de los betas, la medida de desempeño. Esta medida de desempeño es conocida también en la literatura como función de pérdida, queriendo indicar con ello una pérdida de ajuste del modelo a los datos. Es importante aclarar que no existe una única elección para la función de pérdida de un modelo predictivo, otras alternativas típicas son la mediana de las desviaciones absolutas (MEDA), la media de los errores absolutos (MAE), entre otros. Sin embargo en regresión lineal es usual usar MSE porque su expresión como cuadrado de errores permite usar teoría de inferencia basada en la distribución F cuando los errores se asumen normales. Veremos esto más a detalle cuando discutamos pruebas de hipótesis sobre el modelo de regresión.

1.10 Una relación importante entre los betas del modelo lineal

Retomando nuestro modelo original

\[y = \beta_{0} + \beta_{1} * x + e \text{ con } e \text{~}N(0, {\sigma_e}^2)\]

Tomando valor esperado en ambos lados de la ecuación tenemos que:

\[E(y) = E(\beta_{0} + \beta_{1} * x + e)\] Pero como \(\beta_0\) y \(\beta_1\) son constantes, tenemos que:

\[E(y) = E(\beta_0) + E(\beta_1*x) + E(e)\] Así que:

\[E(y) = \beta_0 +\int_{-\infty}^{\infty}\beta_1 *xf(x)dx + E(e)\] \[E(y) = \beta_0 + \beta_1*\int_{-\infty}^{\infty}xf(x)dx + E(e)\] Y como \(\int_{-\infty}^{\infty}xf(x)dx\) es por definición \(E(x)\) y además \(E(e) = 0\) por diseño (de acuerdo con el supuesto usado para la distribución del error), entonces:

\[E(y) = \beta_0 + \beta_1*E(x)\] Dicho de otra manera, la recta del modelo debe pasar por el punto \((E(x), E(y))\) !!. Este punto no es más que el centroide de distribución bivariada \((x_i, y_i)\) de mi proceso generador de datos.

Ciertamente no conocemos el valor esperado ni de la variable predictora X ni de la variable respuesta Y, pero buenos estimados basados en información muestral serían \(\bar{x}\) y \(\bar{y}\). Es decir, no podremos saber con precisión los valores \(\beta_0\) y \(\beta_1\), pero nos podemos conformar con buenas estimaciones que cumplan la condición:

\[\bar{y} = \color{red}{\widetilde{\beta}_0} + \color{red}{\widetilde{\beta}_1}\bar{x}\]

El coloreado y la tilde ancha para los betas pretenden enfatizar que estos valores son estimados de los correspondientes \(\beta_0\) Y \(\beta_1\) del modelo teórico subyacente. En conclusión: la mejor recta estimada deberá pasar por \((\bar{x}, \bar{y})\), y por lo tanto hay una relación atando a los betas estimados:

\[\begin{equation} \label{eq:ec0} \color{red}{\widetilde{\beta}_0} = \bar{y} - \color{red}{\widetilde{\beta}_1}\bar{x} \end{equation} \]

Es decir que dando valores a \(\color{red}{\widetilde{\beta}_1}\) podremos obtener valores para \(\color{red}{\widetilde{\beta}_0}\) y en últimas dibujar diversas rectas en el plano, con el fin de establecer cuál de ellas minimiza el \(MSE_{\text{NO ingenuo}}\).

1.11 Una animación que ilustra la estimación de betas

A continuación presentamos un código en R que permite comprender cómo funcionan numéricamente los métodos de estimación de betas para un modelo lineal simple, específicamente cuando se usa como método de estimación el criterio de minimización del promedio de los cuadrados de los errores:

# Extraemos las variables de interés
x = hass_dolar$precio_dolar
y = hass_dolar$precio_aguacate_kg

# Calculamos medias muestrales para las variables X y Y

x_barra <- mean(x)
y_barra <- mean(y)

# Definimos posibles valores para b1, y por consiguiente para b0
b1_values <- seq(-1, 3, by = 0.08)
b0_values <- y_barra - b1_values * x_barra

# Fijamos límites en X y Y para cada gráfico
xmin = min(x) - 50
xmax = max(x) + 50
ymin = min(y) - 50
ymax = max(y) + 50

# Construimos una función para calcular el MSE en función de los n parejas (x_i, y_i), b1 y b0
MSE <- function(x, y, b0, b1) {
  errors <- y - (b0 + b1*x)
  squared_errors <- errors^2
  mse <- mean(squared_errors)
  return(mse)
}

# Crear un dataframe para almacenar los resultados
list_mse <- vector()

# Calcular el MSE para diferentes valores de B1
for (i in 1:length(b1_values)) {
  mse <- MSE(x, y, b0_values[i], b1_values[i])
  list_mse[i] <- mse
}

mse_df = data.frame(
  B0 = b0_values,
  B1 = b1_values,
  mse = list_mse,
  color_point = 'green'
)

# Capturando el renglón donde ocurre el mínimo mse
pos_min <- which(mse_df$mse==min(mse_df$mse))
b0_aprox <- round(mse_df[pos_min, "B0"],0)
b1_aprox <- mse_df[pos_min, "B1"]

# Impriendo las filas cercanas al mínimo
mse_df$color_point[(pos_min - 4):(pos_min + 4)] <- rep('red',9)
ftable(mse_df[(pos_min - 8):(pos_min + 8), ])

¿Qué hemos logrado hasta acá?

• Demostrar que el valor de \(MSE_{\text{NO ingenuo}}\) depende de los betas elegidos para el modelo.
• Demostrar cómo el valor de \(\color{red}{\widetilde{\beta}_0}\) está atado con el de \(\color{red}{\widetilde{\beta}_1}\) si de verdad queremos que la recta estimada respete la condición de que \(E(e)=0\).
• Generar un dataframe que calcula los valores de \(MSE_{\text{NO ingenuo}}\) para cada par de valores elegidos para \(\color{red}{\widetilde{\beta}_0}\) y \(\color{red}{\widetilde{\beta}_1}\), dando valores a \(\color{red}{\widetilde{\beta}_1}\)
• Estimamos que el \(\color{red}{\widetilde{\beta}_1}\) óptimo está cerca del valor 0.92, y que el \(\color{red}{\widetilde{\beta}_0}\) óptimo está cerca del valor 1513.

Ahora graficaremos una curva que muestre esta relación:

if (!file.exists("beta_vs_mse.gif")){
  # Crear la animación
colors <- c("green", "red")
animation <- ggplot(mse_df, 
                    aes(x = B1, y = mse)) +
  geom_line(aes(x = B1, y = mse,
                color = factor(color_point),
                group = cumsum(c(0, diff(factor(color_point) != 0)))
  )) +
  geom_point(aes(color = factor(color_point),
                 group = cumsum(c(0, diff(factor(color_point) != 0))))) +
  labs(title = "Valores de MSE para cada posible B1, respetando E(e) = 0") +
  geom_text(data = mse_df, 
            aes(label = paste("MSE: ", round(mse,0))), x= 1, y = 2e6,
            color='black'
  ) +
  geom_text(data = mse_df, 
            aes(label = paste("MSE mínimo: ", round(min(mse),0))),
            x= 1, y = min(mse_df$mse) - 1e5, color='red') +
  scale_color_manual(values = colors,
                     labels = c("Subóptima", "Cercana al óptimo"),
                     name = "Soluciones") +
  theme_minimal() +
  transition_reveal(B1)

  animate(animation, duration = 10, fps = 2, renderer = gifski_renderer())
  anim_save("beta_vs_mse.gif")
}

Ahora veamos el efecto que tiene la elección de los betas en la recta regresora, e incluyamos pausas en la animación para detenernos en el momento en que encontremos la recta óptima.

if (!file.exists("beta_estimation.gif")){
  
# Las pausas en la animación son fotogramas con información repetida correspondiente
# al renglón óptimo

fotogramas_pausa <- data.frame(
  B0 = rep(b0_values[pos_min], 2*nrow(mse_df)),
  B1 = rep(b1_values[pos_min], 2*nrow(mse_df)),
  mse = rep(list_mse[pos_min], 2*nrow(mse_df)),
  color_point = rep('red',2*nrow(mse_df))  
)

mse_df <- rbind(mse_df, fotogramas_pausa)
mse_df <- mse_df[order(mse_df$B1), ]
mse_df$estado <- seq(1,nrow(mse_df))

# Construimos el dataframe de predicciones para cada estado, con el cual poder
# graficar los residuales en cada fotograma

df_predicciones <- data.frame(B0_pred=numeric(0),
                              B1_pred=numeric(0),
                              x_from_pred=numeric(0),
                              y_from_pred=numeric(0),
                              y_pred_mod=numeric(0),
                              estado = numeric(0))
for (i in 1:nrow(mse_df)){
  df_nuevo <-data.frame(
    BO_pred = rep(mse_df[i, 'B0'], length(x)),
    B1_pred = rep(mse_df[i, 'B1'], length(x)),
    x_from_pred = x,
    y_from_pred = y,
    y_pred_mod = mse_df[i, 'B0'] + mse_df[i,'B1']*x,
    estado = rep(mse_df[i, 'estado'], length(x))
  ) 
  df_predicciones <- rbind(df_predicciones, df_nuevo)
}

# Creamos el gráfico base

base_plot <- ggplot(data.frame(x = x, y = y), aes(x = x, y = y)) +
  geom_point() +
  geom_point(x=mean(x), y=mean(y), color='green', size=2.5) +
  geom_abline(aes(slope = 0, intercept = y_barra), 
              linetype = "dashed", color = "blue") +
  labs(title = "Estimación de Betas por Minimización del MSE") +
  theme_minimal()

# Creamos la animación

animation <- base_plot +
  geom_abline(aes(slope = B1, intercept = B0, color = color_point), data=mse_df) +
  geom_text(aes(label = paste("MSE: ", round(mse,0))), x= 2000, y = 6000, data=mse_df,
            color = "red")+
  geom_segment(data = df_predicciones,
               aes(x=x_from_pred, y=y_from_pred, xend=x_from_pred, yend=y_pred_mod),
               col = "violet", lty='dashed') +
  transition_states(estado, transition_length = 2, state_length = 1) +
  scale_color_manual(values = c("red" = "red", "green" = "green"),
                     labels = c("Subóptima", "Óptima"),
                     name = "Rectas regresoras") +
  enter_fade() +
  exit_fade()

# Guardamos la animación en un archivo GIF
animate(animation, duration = 20, fps = 2, renderer = gifski_renderer())
anim_save("beta_estimation.gif", animation)
}

1.12 Usando funciones en R para estimar parámetros

Finalmente hemos llegado al punto en que usaremos las funciones de R para estimar parámetros del modelo de regresión lineal. Esto se hará bajo el enfoque de optimización, para ello conviene darle una mirada al siguiente problema:

Estime un par de valores \((\widetilde{\beta}_0,\widetilde{\beta}_1)\) tal que minimice la expresión \(MSE = f(\beta_0,\beta_1)\), donde:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i - \beta_0 - \beta_1*x_i)²\] Tomando en cuenta que no hay restricciones aplicadas a los betas.

Para resolver el problema primero efectuemos algunas simplificaciones partiendo de la siguiente identidad algebraica: \((a+b+c)² = a²+b²+c²+2ab+2ac+2bc\), de tal manera que:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}(y_i²+\beta_0²+\beta_1²x_i²-2y_i\beta_0-2y_i\beta_1x_i+2\beta_0\beta_1x_i)\]

Ahora distribuyamos la suma y separemos los términos que dependen de su índice:

\[f(\beta_0,\beta_1) = \frac{1}{n}\left[\sum^{n}_{i=1}y_i²+\beta_0²\sum^{n}_{i=1}1+\beta_1²\sum^{n}_{i=1}x_i²-2\beta_0\sum^{n}_{i=1}y_i-2\beta_1\sum^{n}_{i=1}y_ix_i+2\beta_0\beta_1\sum^{n}_{i=1}x_i\right]\] Simplificando términos y distribuyendo la fracción:

\[f(\beta_0,\beta_1) = \frac{1}{n}\sum^{n}_{i=1}y_i²+\frac{\beta_0²}{n}n+\frac{\beta_1²}{n}\sum^{n}_{i=1}x_i²-\frac{2\beta_0}{n}\sum^{n}_{i=1}y_i-\frac{2\beta_1}{n}\sum^{n}_{i=1}y_ix_i+\frac{2\beta_0\beta_1}{n}\sum^{n}_{i=1}x_i\] Para finalmente obtener:

\[f(\beta_0,\beta_1) = \beta_0² +\bar{x²}\beta_1² -2\bar{y}\beta_0-\left[\frac{2}{n}\sum^{n}_{i=1}x_iy_i\right]\beta_1+2\bar{x}\beta_0\beta_1+\bar{y²} \] En la expresión anterior, las cantidades \(\bar{x²}\), \(-2\bar{y}\), \(-\frac{2}{n}\sum^{n}_{i=1}x_iy_i\), \(2\bar{x}\), \(\bar{y²}\) son todas constantes por depender sólo de los datos que alimentan al modelo, en ese sentido son coeficientes numéricos para los respectivos betas de la función.

LLegados a este punto, podemos sacar las siguientes conclusiones:

• La función \(f\) sólo depende de los betas, puesto que los demás valores son constantes, es decir los valores \((x_i, y_i)\) son números consignados en una base de datos, y \(n\) es la cantidad de pares de números consignados. Todos los coeficientes dependerán de operaciones sobre estos números.

• De acuerdo con lo anterior la función \(f\) es una función multivariada, pues depende de dos variables: \(\beta_0\) y \(\beta_1\).

• Existen métodos muy sencillos para resolver un problema de optimización no restringido cuando la función a optimizar es convexa, entonces sólo nos resta demostrar que la función anterior lo es.

Definición de convexidad:

Para demostrar que \(f(\beta_0, \beta_1)\) es convexa, debemos verificar que su matriz hessiana sea semidefinida positiva. La matriz hessiana es una matriz de segundas derivadas parciales de la función.

Primero, calculemos las derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):

\[\begin{equation} \label{eq:ec1} \frac{\partial f}{\partial \beta_0} = 2\beta_0-2\bar{y}+2\bar{x}\beta_1 \end{equation} \] \[\begin{equation} \label{eq:ec2} \frac{\partial f}{\partial \beta_1} = 2\bar{x²}\beta_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\beta_0 \end{equation} \]

Luego, calculemos las segundas derivadas parciales de \(f\) con respecto a \(\beta_0\) y \(\beta_1\):

\[\frac{\partial^2 f}{\partial \beta_0^2} = 2 > 0\] \[\frac{\partial^2 f}{\partial \beta_1^2} = 2\bar{x^2} > 0\] \[\frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} = 2\bar{x}\]

La matriz hessiana de \(f\) es:

\[ H(f) = \begin{bmatrix} \frac{\partial^2 f}{\partial \beta_0^2} & \frac{\partial^2 f}{\partial \beta_0 \partial \beta_1} \\ \frac{\partial^2 f}{\partial \beta_1 \partial \beta_0} & \frac{\partial^2 f}{\partial \beta_1^2} \end{bmatrix} \]

Por lo tanto:

\[ H(f) = \begin{bmatrix} 2 & 2\bar{x}\\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \]

Para que \(f(\beta_0, \beta_1)\) sea convexa, la matriz hessiana debe ser semidefinida positiva, lo que significa que debe cumplir la propiedad:

\[v^ \top Hv>=0 \text{ para }\forall v \in \mathbb{R}² \] Es decir, debemos verificar si:

\[\begin{bmatrix} a & b \end{bmatrix}H\begin{bmatrix} a \\ b \\ \end{bmatrix}>=0 \text{ para } \forall (a,b) \]

Veamos si esto se cumple:

\[ \begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} 2 & 2\bar{x} \\ 2\bar{x} & 2\bar{x^2} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\= \begin{bmatrix} 2a+2b\bar{x} & 2a\bar{x} + 2b\bar{x²} \end{bmatrix} \begin{bmatrix} a \\ b \\ \end{bmatrix}\\=2a²+2ab\bar{x}+2ab\bar{x}+2b²\bar{x²} \\=2a²+4ab\bar{x}+2b²\bar{x²}\\ =2(a²+2ab\bar{x}+b²(\bar{x})²)-2b^2(\bar{x})²+2b²\bar{x²}\\ =2(a+b\bar{x})²+2b²(\bar{x²}-(\bar{x})²) \]

Acá es interesante notar que:

\[\bar{x²}-(\bar{x})²=\frac{1}{n}\left[\sum^{n}_{i=1}x_i²\right]-\bar{x}\bar{x} =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}(n\bar{x})\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}x_i²-\bar{x}\sum^{n}_{i=1}x_i\right] =\frac{1}{n}\left[\sum^{n}_{i=1}(x_i²-\bar{x}x_i)\right] \\ =\frac{1}{n}\left[\sum^{n}_{i=1}((x_i²-2\bar{x}x_i)+\bar{x}x_i)\right] =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i²-2\bar{x}x_i+\bar{x}²)+(\bar{x}x_i-\bar{x}²)\right)\right]\\ =\frac{1}{n}\left[\sum^{n}_{i=1}\left((x_i-\bar{x})²+\bar{x}(x_i-\bar{x})\right)\right] =\frac{1}{n}\sum^{n}_{i=1}(x_i-\bar{x})²+\bar{x}\left[\sum^{n}_{i=1}\frac{x_i}{n}-\sum^{n}_{i=1}\frac{\bar{x}}{n}\right]\\ =\frac{S_{xx}}{n}+\bar{x}\left(\bar{x}-\bar{x}\sum^{n}_{i=1}\frac{1}{n}\right)\\ =\frac{S_{xx}}{n}+\bar{x}(\bar{x}-\bar{x}*1)=\frac{S_{xx}}{n}\]

Obteniendo finalmente

\[\begin{equation} \label{eq:sigma_x} \bar{x²}-(\bar{x})² = \frac{S_{xx}}{n} \end{equation} \] Note que aquí hemos usado la siguiente definición:

\[S_{xx}=\sum^{n}_{i=1}(x_i-\bar{x})*(x_i-\bar{x})=\sum^{n}_{i=1}(x_i-\bar{x})^2\] Luego:

\[v^ \top Hv = 2(a+b\bar{x})²+2b²\frac{S_{xx}}{n} \text{ con } \text{ }v=(a,b)\] Y como ambos sumandos son no negativos se tiene que \[v^ \top Hv >= 0 \text{ para }\forall v\in \mathbb{R}²\] Es decir, la matriz hessiana es semidefinida positiva y por lo tanto \(f(\beta_0,\beta_1)\) es convexa en el espacio de parámetros \(\beta_0\) y \(\beta_1\). Esto implica que el problema de regresión lineal, que busca minimizar esta función de errores, tiene un mínimo global y puede ser resuelto de manera eficiente mediante métodos de optimización convexa.

Una condición necesaria y suficiente para que \((\widetilde{\beta}_0,\widetilde{\beta}_1)\) sea un mínimo global es que \(\nabla{} f(\widetilde{\beta}_0,\widetilde{\beta}_1)=0\)

Por lo tanto retomando \(\ref{eq:ec1}\) y \(\ref{eq:ec2}\) el problema de minimización queda resuelto hallando la solución al sistema de ecuaciones:

\[\begin{equation} \label{eq:ec3} 2\widetilde{\beta}_0-2\bar{y}+2\bar{x}\widetilde{\beta}_1=0 \end{equation} \] \[\begin{equation} \label{eq:ec4} 2\bar{x²}\widetilde{\beta}_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}\widetilde{\beta}_0=0 \end{equation} \]

De modo que despejando \(\widetilde{\beta}_0\) en \(\ref{eq:ec3}\) tenemos lo siguiente:

\[\begin{equation} \label{eq:ec5} \widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x} \end{equation} \] Esto en concordancia con \(\ref{eq:ec0}\), y además haciendo reemplazos en \(\ref{eq:ec4}\) obtenemos:

\[ 2\bar{x²}\widetilde{\beta}_1-\frac{2}{n}\sum^{n}_{i=1}x_iy_i+2\bar{x}(\bar{y}-\widetilde{\beta}_1\bar{x})=0\\ 2\bar{x²}\widetilde{\beta}_1-2\bar{x}²\widetilde{\beta}_1 = \frac{2}{n}\sum^{n}_{i=1}x_iy_i-2\bar{x}\bar{y}\\ (\bar{x²} - \bar{x}²)\widetilde{\beta}_1 = \frac{1}{n}\sum^{n}_{i=1}x_iy_i-\bar{x}\bar{y}\\\]

Y usando la expresión dada en \(\ref{eq:sigma_x}\) tenemos:

\[\frac{S_{xx}}{n}\widetilde{\beta}_1=\frac{1}{n}\left(\sum^{n}_{i=1}x_iy_i-n\bar{x}\bar{y}\right)\\S_{xx}\widetilde{\beta}_1 = \sum^{n}_{i=1}x_iy_i-n\bar{x}\bar{y} \]

Ahora, podemos usar el hecho de que el lado derecho de la ecuación anterior cumple la siguiente equivalencia:

\[\sum^{n}_{i=1}x_iy_i-n\bar{x}\bar{y}= \sum^{n}_{i=1}(x_i-\bar{x})*(y_i-\bar{y}) = S_{xy}\] Es decir, se puede convertir en la suma de productos cruzados de la variable \(x\) y la variable respuesta \(y\), previamente centradas respecto a sus medias.

Por lo tanto se tendría que:

\[\begin{equation} \label{eq:ec6} \widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}} \end{equation} \]

Las ecuaciones \(\ref{eq:ec5}\) y \(\ref{eq:ec6}\) son muy informativas en su estructura, la primera indica que la recta regresora debe pasar por el centroide de la nube de puntos: \((\bar{x}, \bar{y})\) y la segunda indica que la pendiente de la recta es proporcional al grado de relación lineal que tenga la variable \(x\) con la variable \(y\) medida a través de la suma de productos cruzados \(S_{xy}\), la cual a su vez se relaciona con la covarianza entre las variables \(x\) y \(y\).

Estos sencillos cálculos son realizados por R a través de la función “lm”, así:

# b1 = cov(x,y) / (sigma_x)²
b1 <- cov(x,y) / var(x)
print(b1)

## [1] 0.9449775

# bo = y_barra - b1*x_barra
b0 <- mean(y) - b1*mean(x)
print(b0)

## [1] 1436.679

residuales <- y - (b0 + b1*x)
n <- length(x)

sigma_2 <- sum(residuales^2)/(n-2)
print(sigma_2)

## [1] 288485.9

# Confirmamos usando la función lm, a continuación la documentación

# lm(formula, data, subset, weights, na.action,
#    method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
#    singular.ok = TRUE, contrasts = NULL, offset, ...)

mod1 <- hass_dolar %>% lm(precio_aguacate_kg ~ precio_dolar,.)
print(mod1$coefficients)

##  (Intercept) precio_dolar 
## 1436.6790256    0.9449775

resumen <-  summary(mod1)
resumen$sigma

## [1] 537.1089

parametros <- data.frame(metodo = c("manual","lm"), 
                         beta_0 = c(b0,mod1$coefficients[1]),
                         beta_1 = c(b1,mod1$coefficients[2]),
                         sigma_2 = c(sigma_2, (resumen$sigma)^2))



parametros %>% 
  ftable()

#print(mod1$)
# Generando el modelo con R Base sería así:
# mod1 <- lm(promedio_aguacate_kg ~ promedio_dolar, hass_dolar)

Podemos notar entonces que ambos métodos son consistentes entre si, algo que no es de extrañar pues internamente R implementa estos métodos en sus librerías.

6.1 Estudiando el comportamiento estadístico de los estimadores

¿Cómo garantizamos que dentro de todas las posibles operaciones matemáticas que podemos ejecutar con los valores de la muestra, las operaciones:

\[\widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}}\] \[\widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x}\] Sean las operaciones que conducen a valores más próximos a \(\beta_1\) y \(\beta_0\) respectivamente, aun si la muestra de datos cambia?.

No podemos demostrar esto con el rigor necesario porque escapa un poco del alcance del presente curso, pero podemos al menos ilustrar lo que queremos indicar por cercanía a los betas del modelo usando el siguiente enfoque analítico.

Detengamonos un instante a analizar la estructura de \(\widetilde{\beta}_1\), notemos que este estimador puede pensarse como una combinación lineal de las \(y_i\), veamos por que:

\[\widetilde{\beta}_1 = \frac{S_{xy}}{S_{xx}}=\frac{\sum_{i=1}^{n}(x_i-\bar{x})*(y_i-\bar{y})}{S_{xx}}=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})-\bar{y}(x_i-\bar{x})}{S_{xx}}\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\sum_{i=1}^{n}\frac{\bar{y}(x_i-\bar{x})}{S_{xx}}\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\sum_{i=1}^{n}(x_i-\bar{x})\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\left(\sum_{i=1}^{n}x_i-\sum_{i=1}^{n}\bar{x}\right)\\=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}-\frac{\bar{y}}{S_{xx}}\left(n\bar{x}-n\bar{x}\right)=\sum_{i=1}^{n}\frac{y_i(x_i-\bar{x})}{S_{xx}}\] Ahora sea \(c_i= \frac{x_i-\bar{x}}{S_{xx}}\) entonces concluimos que:

\[\widetilde{\beta}_1=\sum_{i=1}^{n}c_i*y_i\] Es decir \(\widetilde{\beta}_1\) es en efecto una combinación lineal de los \(y_i\) con coeficientes \(c_i\), por lo tanto, retomando a nuestro proceso generador de datos, y tomando un \(y_i\) particular generado por dicho proceso, tenemos que \(y_i = \beta_0 + \beta_1x_i + e_i\) donde los \(e_i\) se asumen normales e independientes, lo cual implica que:

\[E(y_i|x_i) = E((\beta_0 + \beta_1x_i + e_i)|x_i)=\beta_0+\beta_1x_i + E(e_i) = \beta_0 + \beta_1x_i\]

Ya que \(E(e_i)=0\) desde los supuestos. Y por otro lado que:

\[Var(y_i|x_i) = Var((\beta_0 + \beta_1x_i + e_i)|x_i) = 0+0+Var(e_i)= {\sigma_e}^2\]

Ya que \(Var(e_i) = {\sigma_e}^2\) también desde los supuestos.

Así que como conclusión tenemos que, cuando el \(x_i\) está ya dado:

\[y_i\text{~}N(\beta_0 + \beta_1x_i,{\sigma_e}^2)\]

Todo esto por propiedad de linealidad de la normalidad, es decir, debido a que una transformación afín de la variable normal \(e_i\) dada por un desplazmiento \(\beta_0\) y un factor de escala \(\beta_1\) es también una variable aleatoria normal.

Lo anterior, tácitamente nos está diciendo que nuestro estimador \(\widetilde{\beta}_1\) es entonces una combinación lineal de variables aleatorias normales, razón por la cual, partiendo de la propiedad reproductiva de la normal podremos afirmar que, cuando todos los \(x_i\) están dados \(\widetilde{\beta}_1\) es también una variable aleatoria normal.

Investiguemos su media:

\[E(\widetilde{\beta}_1)=E\left(\sum_{i=1}^{n}c_iy_i\right)=\sum_{i=1}^{n}E(c_iy_i)=\sum_{i=1}^{n}c_iE(y_i)=\sum_{i=1}^{n}c_i(\beta_0 + \beta_1x_i)\\=\beta_0\sum_{i=1}^{n}c_i + \beta_1\sum_{i=1}^{n}c_ix_i\] Pero es fácil ver que \(\sum_{i=1}^{n}c_i=0\) y que \(\sum_{i=1}^{n}c_ix_i=1\) (por favor demuéstrelo), y entonces:

\[E(\widetilde{\beta}_1)=\beta_1\]

Lo que indicaría que \(\widetilde{\beta}_1\) es un estimador insesgado para \(\beta_1\)!!!

Y ahora investiguemos su varianza:

\[Var(\widetilde{\beta}_1)=Var\left(\sum_{i=1}^{n}c_iy_i\right)=\sum_{i=1}^{n}{c_i}²Var(y_i)=\sum_{i=1}^{n}{c_i}²{\sigma_e}²={\sigma_e}²\sum_{i=1}^{n}{c_i}²\] Ya que de nuevo por asunciones sobre el proceso generador de datos las \(y_i\) son no correlacionadas y habíamos demostrado previamente que \(Var(y_i) = {\sigma_e}²\), sin importar el \(y_i\) considerado, es decir, estos son homocedásticos como consecuencia de la definición del modelo. Por lo tanto:

\[Var(\widetilde{\beta}_1)={\sigma_e}²\sum_{i=1}^{n}{\left(\frac{x_i-\bar{x}}{S_{xx}}\right)}^2=\frac{\sigma_e²}{{S²}_{xx}}\sum_{i=1}^{n}(x_i-\bar{x})²=\frac{\sigma_e²}{{S²}_{xx}}*S_{xx}=\frac{\sigma_e²}{{S}_{xx}}\] Finalmente, ya que hemos investigado ampliamente el comportamiento del estimador \(\widetilde{\beta}_1\) usemos este conocimiento para investigar el comportamiento de \(\widetilde{\beta}_0\). Para ello, considere que \(\widetilde{\beta}_0 = \bar{y} - \widetilde{\beta}_1\bar{x}\), por lo tanto:

\[E(\widetilde{\beta}_0)=E(\bar{y} - \widetilde{\beta}_1\bar{x})=E(\bar{y})-\bar{x}E(\widetilde{\beta}_1)=E(\sum_{i=1}^{n}\frac{y_i}{n})-\bar{x}\beta_1=\frac{1}{n}\sum_{i=1}^{n}E(y_i)-\bar{x}\beta_1\\=\frac{1}{n}\sum_{i=1}^{n}(\beta_0 + \beta_1x_i)-\bar{x}\beta_1=\frac{1}{n}\left(\sum_{i=1}^{n}\beta_0+\beta_1\sum_{i=1}^{n}x_i\right)-\bar{x}\beta_1\\=\frac{1}{n}(n\beta_0 + \beta_1n\bar{x})-\bar{x}\beta_1=\beta_0\]

De nuevo \(\widetilde{\beta}_0\) es un estimador insesgado de \(\beta_0\)!!!

Para la varianza tenemos:

\[Var(\widetilde{\beta}_0)=Var(\bar{y} - \widetilde{\beta}_1\bar{x})=Var(\bar{y})+{\bar{x}}^2Var(\widetilde{\beta}_1)=\frac{{\sigma_e}²}{n} + {\bar{x}}^2\frac{{\sigma_e}²}{S_{xx}}={\sigma_e}²\left(\frac{1}{n}+\frac{{\bar{x}}^2}{S_{xx}}\right)\]

También es importante hacer ver que \(\widetilde{\beta}_0\) puede pensarse como combinación lineal de las \(y_i\) ya que en su fórmula están involucradas las cantidades \(\bar{y}\) y \(\widetilde{\beta}_1\) a través de una combinación lineal. Por lo tanto, como una combinación lineal de combinaciones lineales de las \(y_i\) es también una combinación lineal de \(y_i\) podemos afirmar que por propiedad reproductiva de la normal \(\widetilde{\beta}_0\) también se distribuye normal.

Una última anotación a realizar es que según el teorema de Gauss-Markov, para el modelo de regresión lineal con los supuestos \(E(e)=0\) y \(Var(e)={\sigma_e}²\) y por supuesto con errores no correlacionados, los estimadores por mínimos cuadrados son insesgados (tal como acabamos de demostrar), pero más importante que eso, tienen la varianza más pequeña que se pueda esperar dentro del universo de estimadores insesgados que se construyen como combinaciones lineales de los \(y_i\). Es por esta razón que con frecuencia se dice que los estimadores por mínimos cuadrados son los estimadores lineales insesgados óptimos, donde óptimos implica que son de varianza mínima.

6.2 Construyendo estadísticos de prueba para hipótesis sobre los beta

Tanto para \(\widetilde{\beta}_0\) como para \(\widetilde{\beta}_1\) hemos deducido expresiones para su varianza y su media, y además hemos demostrado analíticamente que ambos estimadores son variables aleatorias normales. La única limitación que tenemos hasta ahora para realizar pruebas de hipótesis sobre sus valores es el hecho de que la fórmula de la varianza de los betas está involucrando la varianza \({\sigma^2}_e\) la cual es desconocida.

Sin embargo, si inspeccionamos la estructura de la suma de cuadrados de los residuales, tenemos que:

\[SS_{Res} =\sum^{n}_{i=1}(y_i - \hat{y_i})² =\sum^{n}_{i=1}(y_i - \widetilde{\beta}_0 -\widetilde{\beta}_1x_i)²\]
Y como esta suma de cuadrados tiene \(n-2\) grados de libertad, porque dos grados de libertad se asocian con \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\), es posible demostrar que \(E(SS_{Res})=(n-2)*{{\sigma_e}^2}\), de modo que la cantidad:

\[\widetilde{\sigma}²_e = MS_{Res} = \frac{SS_{Res}}{n-2}\] Sería un estimador insesgado para \({\sigma^2}_e\).

También es posible demostrar que la cantidad:

\[\chi²=\frac{SS_{Res}}{{\sigma²}_e}=\frac{(n-2)\widetilde{\sigma}²_e}{{\sigma^2}_e}\]

se distribuye como una variable aleatoria chi-cudrada con \(n-2\) grados de libertad, pues precisamente es la suma de \(n-2\) cuadrados de variables aleatorias normales independientes e idénticamente distribuidas divididas entre su varianza común. También note que este comportamiento teórico esperado para el estimador de \({\sigma²}_e\) depende fuertemente de los residuales del modelo, es decir, puede verse afectado seriamente si el comportamiento de los residuales no se cumple, esto es: independencia, homocedasticidad y normalidad.

Prueba de hipótesis para beta_0:

Considere la siguiente cantidad:

\[z_0 = \frac{\widetilde\beta_0 - \beta_{00}}{\sqrt{Var(\widetilde\beta_0)}}=\frac{\widetilde\beta_0 - \beta_{00}}{\sqrt{{\sigma_e}²\left(\frac{1}{n}+\frac{{\bar{x}}^2}{S_{xx}}\right)}}\] Dada la normalidad de \(\widetilde{\beta_0}\) y el hecho de que \(\widetilde{\beta_0}\) es un estimador insesgado de \(\beta_0\) podemos decir que \(\z_0\) se distribuye como una variable normal estándar. Pero como obviamente tampoco conocemos el verdadero valor de la varianza de \(\widetilde{\beta}_0\), pues está dependiendo de \({\sigma_e}²\) en su lugar tendremos que proponer un estadístico de prueba alternativo dado por:

\[t_0 = \frac{\widetilde\beta_0 - \beta_{00}}{\sqrt{\frac{SS_{Res}}{n-2}\left(\frac{1}{n}+\frac{{\bar{x}}^2}{S_{xx}}\right)}}\] Note como \({\sigma_e}²\) ha sido sustituido por su estimador insesgado deducido previamente. Esta nueva cantidad se distribuirá, según la teoría estadística, como una variable aleatoria t-student con \(n-2\) grados de libertad!!!

Entonces podemos proponer el test de hipótesis:

\(H_0: \beta_0 = \beta_{00}\) \(H_1: \beta_0 \ne \beta_{00}\)

Donde la hipótesis nula se rechaza si \(|t_0|=t_{\frac{\alpha}{2},n-2}\)

Prueba de hipótesis para beta_1:

Similar a como se procedió para \(\beta_0\), para \(\beta_1\) podemos usar el estadístico de prueba dado por:

\[t_0 = \frac{\widetilde\beta_1 - \beta_{10}}{\sqrt{\frac{SS_{Res}}{(n-2){S}_{xx}}}}\]

Y entonces podemos proponer el test de hipótesis:

\(H_0: \beta_1 = \beta_{10}\) \(H_1: \beta_1 \ne \beta_{10}\)

Donde la hipótesis nula se rechaza si \(|t_0|=t_{\frac{\alpha}{2},n-2}\)

Valga aclarar que, como en el caso del estimador \(\widetilde{\sigma}²_e\), el comportanmiento distribucional esperado para los estadísticos de prueba de \(\beta_0\) y \(\beta_1\) puede resultar seriamente comprometido si los residuales no cumplen los supuestos del modelo.

6.3 Simulando el comportamiento estadístico de los estimadores

Ahora podemos comprobar cómo se comportan los valores de \(\widetilde{\beta}_0\), \(\widetilde{\beta}_1\) y \(MS_{Res}\) a lo largo de diversas muestras de tamaño \(n\) obtenidas con el modelo teórico. El pseudo código es:

1. Generar una función para extraer muestras de tamaño \(n\) bajo un modelo teórico dado.
2. Generar una función para ajustar una recta a la muestra anterior usando los estimadores \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) dados en las ecuaciones anteriores.
3. Dibujar algunas rectas ajustadas para un conjunto de 5 muestras obtenidas del modelo teórico con \(\sigma=4\).
4. Dibujar algunas rectas ajustadas para un conjunto de 5 muestras obtenidas del modelo teórico con \(\sigma=10\).
5. Concluir según los resultados.

# Definimos un modelo teórico
n <- 50
vector_x <- c(1:n) * 30/50
sigma <- 4
beta_0 <- 2
beta_1 <- 0.5

# Creamos la función generadora de muestras
generador_muestra <-  function(n,vector_x,sigma, beta_0,beta_1){
  error <- rnorm(n,mean = 0, sd = sigma)
  y_observada <- beta_0 + beta_1 * vector_x + error
  return( list(x = vector_x, y = y_observada) )
}

# Creamos la función que estima el error estándar de b0 y b1
se_betas <- function(vector_x, residuales, n){
  SS_res <-sum(residuales^2)
  MS_res <- SS_res / (n-2)
  x_barra <- mean(vector_x)
  SS_x <- sum((vector_x-x_barra)^2)
  se_b0 <- sqrt(MS_res*((1/n) + (x_barra^2)/SS_x))
  se_b1 <- sqrt(MS_res/SS_x)
  return(list(se_b0 = se_b0, se_b1 = se_b1))
}

# Corroborando que la función quedó bien 👍
generador_muestra(10,c(1:10),sigma, beta_0,beta_1)

## $x
##  [1]  1  2  3  4  5  6  7  8  9 10
## 
## $y
##  [1]  2.7913467 -4.6325215  3.1425105  0.2917908  3.2684051  4.0165954
##  [7]  2.2582931  5.3461087  4.3314844 10.2005543

# Creamos la función que estima una recta sobre una muestra dada
estimador_betas <- function(x,y){
    model <- lm(y~x)
    coeficientes <- model$coefficients
    errores_betas <- se_betas(x, model$residuals, length(x))
    return(list(beta_0 = coeficientes[1],
                beta_1 = coeficientes[2],
                se_b0 = errores_betas$se_b0,
                se_b1 = errores_betas$se_b1,
                sigma = sigma(model)
                ))  
}

# Ajustamos semilla para reproducibilidad y una función que genera un
# gráfico con 5 rectas estimadas

set.seed(123)

generar_grafico <- function(x,n,sigma,beta_0,beta_1,position_legend){
df_rectas <- data.frame(x=numeric(),
                        y_pred=numeric(),
                        numero_muestra = numeric())
df_true <- data.frame(x=x,
                      y_true = (beta_0 + beta_1 * x))
for (i in 1:5) {
  muestra <- generador_muestra(n,x,sigma, beta_0,beta_1)
  betas_estimados <-  estimador_betas(muestra$x, muestra$y)
  df <- data.frame(x = x,
                   y_pred = (betas_estimados$beta_0 + 
                          betas_estimados$beta_1 * x),
                   numero_muestra = rep(i,n))
  df_rectas <- rbind(df_rectas,df)
}

p <- ggplot() +
  geom_line(data = df_rectas, aes(x=x, y=y_pred,
                                  color = factor(numero_muestra)))+
  scale_color_manual(values = c("1" = "blue", "2" = "black",
                                "3" = "green", "4" = "violet",
                                "5" = "orange",
                                "y real" = "red"),
                     name = "Rectas estimadas",
                     labels = c(paste("Muestra_",seq(1,5)),
                                "Recta real"))+
  labs(x="Valor x",
       y="y predicho",
       title = "Rectas estimadas sobre muestras distintas",
       subtitle = paste0("Sigma=",sigma, " n=",n))+
  geom_point(data = df_true, aes(x=x, y=y_true, color="y real"))+
  theme(legend.position=position_legend)
  return(p)
}

# Invocamos la función con valores de sigma distintos (4 y 10) y
# presentamos un gráfico comparativo.

grid.arrange(generar_grafico(x=vector_x,n=50,sigma=4,beta_0,beta_1,
                             position_legend='right'),
             generar_grafico(x=vector_x,n=50,sigma=10,beta_0,beta_1,
                             position_legend='right'),
             nrow=2)

Conclusiones:

• Las rectas estimadas se aproximan a la recta real, pero dicha cercanía dependerá del sigma (variabilidad inherente del proceso generador de datos)
• Como las rectas se estan distanciando de la recta teórica conviene investigar qué lo provoca, esto es, cuánto se distancian los betas estimados de los betas teóricos.
• Conviene también estudiar el comportamiento estadístico de los betas estimados, no sólo la distancia a la que quedan de los betas teóricos.

Observación:

Tal como se aprecia en el gráfico, las rectas se distancian más de la teórica conforme el \(\sigma\) desconocido del proceso crece. Y como vimos en el apartado anterior:

\[\widetilde{\sigma}²_e=MS_{Res} = \sum_{i=1}^{n}\frac{(y_i-\tilde{y_i})²}{n-2}\]

Acá el número 2 que se resta en el denominador representa la cantidad de parámetros en mi modelo,es decir que el denominador debe ser los grados de libertad de la suma de cuadrados de los residuales para garantizar la insesgadez del estimador \(\widetilde{\sigma}²_e\), entendiendo por esto la propiedad que tiene de que su promedio sobre todas las posibles muestras extraídas del proceso generador de datos sea igual a \(\sigma²\), esto es \(E(\widetilde{\sigma}²_e)={\sigma_e}²\)

Ahora estamos en condiciones de investigar el comportamiento estadístico de \(\widetilde{\beta}_0\), \(\widetilde{\beta}_1\), y \(\widetilde{\sigma}²_e\) para ello procederemos bajo el siguiente pseudo-código.

• Generar una lista de tamaños de muestra a recorrer.
• Generar un modelo teórico a usar.
• Recorrer a través de un ciclo for todos los tamaños presentes en la lista.
• Invocar 10000 veces el generador de la muestra sobre cada tamaño de muestra.
• Sobre cada muestra invocar el estimador de los betas y sigma.
• Almacenar cada grupo de estimaciones en sus correspondientes vectores de resultados.
• Conformar un dataframe con el tamaño de muestra y estas estimaciones, al cual añadiremos además una columna con la cantidad dada por \(\frac{(n-2)\widetilde{\sigma}²_e}{{\sigma²}_e}\) la cual se espera distribuida como una chi-cuadrado con \(n-2\) grados de libertad, y con las estandarizaciones de los betas estimados, las cuales se esperan distribuidas como t-student con \(n-2\) grados de libertad.
• Culminado el ciclo producir gráficos de las cantidades de interés: sigma estimado, cantidad ji-cudrado, betas estimados, y betas estandarizados.
• Concluir sobre la validez de los comportamientos distribucionales esperados.
• Concluir sobre cómo el tamaño de muestra afecta la magnitud del error de estimación de los betas.

# Parámetros para fijar el proceso generador de datos
sigma <- 4
beta_0 <- 2
beta_1 <- 0.5

# vectores donde guardaremos los resultados de los betas estimados para
# cada tamaño de muestra
vector_beta_0_est <- vector()
vector_beta_1_est <- vector()
vector_se_beta_0_est <- vector()
vector_se_beta_1_est <- vector()
vector_sigma_est <- vector()
tamanos <- vector()

# Ciclo de alimentación de los vectores

if(!file.exists("df_estimaciones.csv")){
j<-1
for (n in c(5,10,15,20,50,100,500,1000)){
  for (i in 1:10000){
    muestra <-generador_muestra(n=n,vector_x=seq(1,n),sigma=sigma,
                                beta_0=beta_0,beta_1=beta_1)
    betas_estimados <- estimador_betas(muestra$x, muestra$y)
    vector_beta_0_est[j] <- betas_estimados$beta_0
    vector_beta_1_est[j] <- betas_estimados$beta_1
    vector_se_beta_0_est[j] <- betas_estimados$se_b0
    vector_se_beta_1_est[j] <- betas_estimados$se_b1
    vector_sigma_est[j] <- betas_estimados$sigma
    tamanos[j] <- n
    j <- j+1
  }
}

# Construyo el dataframe de estimaciones y en el caso de sigma estimado
# cálculo el estadístico de prueba que se usa para pruebas de hipótesis
# sobre sigma.

df_estimaciones <- data.frame(
  beta_0_est = vector_beta_0_est,
  beta_1_est = vector_beta_1_est,
  beta_0_estandarizado = ((vector_beta_0_est-beta_0) /
                             vector_se_beta_0_est),
  beta_1_estandarizado = ((vector_beta_1_est-beta_1) /
                             vector_se_beta_1_est),
  sigma_est = vector_sigma_est,
  chi_2 = (tamanos -2)*(vector_sigma_est^2)/sigma^2,
  tamano = tamanos
)

write.csv(df_estimaciones,file="df_estimaciones.csv",
          row.names = FALSE, sep=",")
}else{
  df_estimaciones <- read.csv("df_estimaciones.csv", sep=",")
}

df_estimaciones %>% 
  head(10) %>% 
  ftable()

p1 <- ggplot(data=df_estimaciones, aes(x=beta_0_est))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 0 estimado")+
  facet_wrap(vars(tamano), scales="free")

plot(p1)

p2 <- ggplot(data=df_estimaciones, aes(x=beta_0_estandarizado))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 0 estandarizado")+
  facet_wrap(vars(tamano), scales="free")

plot(p2)

p3 <- ggplot(data=df_estimaciones, aes(x=beta_1_est))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 1 estimado")+
  facet_wrap(vars(tamano), scales="free")

plot(p3)

p4 <- ggplot(data=df_estimaciones, aes(x=beta_1_estandarizado))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Beta 1 estandarizado")+
  facet_wrap(vars(tamano), scales="free")

plot(p4)

p5 <- ggplot(data=df_estimaciones, aes(x=sigma_est))+
  geom_histogram()+
  labs(y = "Frecuencia", x="Sigma estimado")+
  facet_wrap(vars(tamano), scales="free")

plot(p5)

p6 <- ggplot(data=df_estimaciones, aes(x=chi_2))+
  geom_histogram()+
  labs(y = "Frecuencia", x="(n-2)^MS_Res/Sigma^2")+
  facet_wrap(vars(tamano), scales="free")

plot(p6)

Conclusiones:

• Los histogramas para \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) tienen forma acampanada lo que sugiere una relación con la distribución t-student. Esto se confirma al inspeccionar los histogramas de las betas estandarizadas los cuales tienden a la normal conforme los grados de libertad crecen. Como una forma de confirmar la distribución subyacente se recomienda al alumno realizar una prueba de kolmogorov sobre cada histograma contra la teórica esperada (una distribución t con n-2 grados de libertad)

• Las varianzas de \(\widetilde{\beta}_0\) y \(\widetilde{\beta}_1\) dependen de los tamaños de muestra (y del sigma del proceso generador de datos que para todos los histogramas se dejó fijo en \(\sigma_e=4\). Puede notar con el rango de variación de los betas se va estrechando conforme el tamaño de muestra crece.

• El cociente \(\frac{(n-2)\widetilde{\sigma}_e^2}{{\sigma_e}^2}\) en efecto parece seguir una distribución chi cuadrado con n-2 grados de libertad. Para confirmarlo se sugiere al alumno aplicar una prueba de kolmogorov sobre cada histograma. Note que por teorema central del límite incluso la ji-cuadrado tenderá a una distribución normal conforme el tamaño de muestra crece.

6.4 Tabla Anova para probar desempeño general del modelo

Finalmente nos resta introducir la tabla de análisis de varianza, o tabla ANOVA. Esta tabla busca realizar un comparativo entre la parte de la varianza de mi variable a predecir \(Y\) que está siendo capturada por mi modelo, vs la parte de la varianza que el modelo no logra capturar. Así pues consideremos la siguiente igualdad:

\[y_i - \bar{y} = (\hat{y}_i - \bar{y}) + (y_i - \hat{y}_i)\] Discutamos cada término:

• \(y_i - \bar{y}\) representa la distancia de cada dato al promedio histórico, es decir mide la variabilidad de la variable original.

• \(\hat{y}_i - \bar{y}\) representa un pedazo de la distancia entre \(y_i\) y \(\bar{y}\), es decir, representa un pedazo de la variabilidad original de \(Y\), ¿Cuál es este pedazo? Aquel que puede ser explicado como consecuencia de la predicción \(\hat{y}_i\). A nosotros nos interesa que este pedazo sea lo más grande que se pueda, precisamente porque queremos que nuestro modelo logre explicar una buena porción de las variaciones que están sufriendo las \(y_i\).

• \(y_i - \hat{y}_i\) es el residual, representa el otro pedazo de la distancia entre \(y_i\) y \(\bar{y}\), es decir, el resto de la variabilidad de \(Y\) que no logró ser capturada por el modelo y que por lo tanto constituye un error de pronóstico. Queremos que esta cantidad sea lo más pequeña posible.

Ahora bien, cómo estas cantidades vienen con signo dependiendo de si la recta predicha se ubica por debajo o por encima de la recta horizontal ingenua, lo mejor es evitar la compensación de errores considerando los cuadrados, así:

\[(y_i - \bar{y})^2 = \left((\hat{y}_i - \bar{y}) + (y_i - \hat{y}_i)\right)²\] Y resolviendo el cuadrado de la derecha tenemos:

\[(y_i - \bar{y})^2 = (\hat{y}_i - \bar{y})^2 + 2(\hat{y}_i - \bar{y})(y_i - \hat{y}_i)+(y_i - \hat{y}_i)^2\] Luego, considerar que en realidad tenemos \(n\) cantidades involucradas acá, una por cada posible \(i=1,..,n\), apliquemos sumatoria en ambos lados para obtener:

\[\begin{equation} \label{eq:sct} \sum^{n}_{i=1}(y_i - \bar{y})^2 = \sum^{n}_{i=1}(\hat{y}_i - \bar{y})^2 + 2\sum^{n}_{i=1}(\hat{y}_i - \bar{y})(y_i - \hat{y}_i)+\sum^{n}_{i=1}(y_i - \hat{y}_i)^2 \end{equation}\] Y en la ecuación \(\ref{eq:sct}\) nos interesa especialmente el término que contiene el doble producto, veamos…

\[2\sum^{n}_{i=1}(\hat{y}_i - \bar{y})(y_i - \hat{y}_i)=2\sum^{n}_{i=1}(\hat{y}_i- \bar{y})e_i=2\sum^{n}_{i=1}\hat{y}_ie_i-2\bar{y}\sum^{n}_{i=1}e_i\] Pero sabemos que \(\sum^{n}_{i=1}e_i=0\) pues los errores se compensan entre sí. Investiguemos el término restante:

\[\sum^{n}_{i=1}\hat{y}_ie_i=\sum^{n}_{i=1}(\widetilde{\beta}_0+\widetilde{\beta}_1x_i)(y_i-(\widetilde{\beta}_0+\widetilde{\beta}_1x_i))\\=\sum^{n}_{i=1}(\bar{y}-\widetilde{\beta}_1\bar{x}+\widetilde{\beta}_1x_i)(y_i-(\bar{y}-\widetilde{\beta}_1\bar{x}+\widetilde{\beta}_1x_i))\\=\sum^{n}_{i=1}[\bar{y}+\widetilde{\beta}_1(x_i-\bar{x})][(y_i-\bar{y})-\widetilde{\beta}_1(x_i-\bar{x})]\\=\bar{y}\sum^{n}_{i=1}(y_i-\bar{y})-\widetilde{\beta}_1\bar{y}\sum^{n}_{i=1}(x_i-\bar{x})+\widetilde{\beta}_1\sum^{n}_{i=1}(x_i-\bar{x})(y_i-\bar{y})-{\widetilde{\beta}_1}^2\sum^{n}_{i=1}(x_i-\bar{x})²\]

Acá es fácil ver que los dos primeros sumandos son nulos pues representan las sumas de las desviaciones de cada variable respecto a su media, luego:

\[\sum^{n}_{i=1}\hat{y}_ie_i=\frac{S_{xy}}{S_{xx}}*S_{xy} - \left(\frac{S_{xy}}{S_{xx}}\right)^2*S_{xx}=0\] Así que retomando de nuevo \(\ref{eq:sct}\) concluimos que:

\[\sum^{n}_{i=1}(y_i - \bar{y})^2 = \sum^{n}_{i=1}(\hat{y}_i - \bar{y})^2 +\sum^{n}_{i=1}(y_i - \hat{y}_i)^2\\ SS_T = SS_R+SS_{Res}\]

Es decir, que en definitiva la varibilidad de \(Y\) capturada a través de \(SS_T\) puede ser descompuesta como la variabilidad explicada por el modelo de regresión (en adelante denominada \(SS_R\)) más la variabilidad de los residuales(en adelante denominada \(SS_{Res}\)). También se ha dicho que \(SS_{Res}\) es tal que \(\frac{SS_{Res}}{{\sigma^2}_e}\) se distribuye como una \({\chi^2}_{n-2}\). Y de manera similar se puede demostrar que \(\frac{SS_{R}}{{\sigma^2}_e}\) se distribuye como una \({\chi^2}_{1}\). Ahora estableciendo sin prueba (ver bibliografía de Montgomery apéndice C3) que ambas sumas de cuadrados son independientes entre sí, se puede concluir que la cantidad:

\[F_0=\frac{\frac{SS_R}{1*{\sigma^2}_e}}{\frac{SS_{Res}}{(n-2)*{\sigma^2}_e}} = \frac{\frac{{\chi^2}_{1}}{1}}{\frac{{\chi^2}_{n-2}}{n-2}}=\frac{MS_{R}}{MS_{Res}}\] Se distribuye como una \(F_{1, n-2}\) de Snedecor con \(1\) grado de libertad en el numerador y \(n-2\) en el denominador, debido a que es el cociente de dos variables aleatorias ji-cuadrado independientes divididas entre sus grados de libertad. Note que hemos introducido una nueva notación para los cocientes \(\frac{SS_R}{1}\) y \(\frac{SS_{Res}}{n-2}\), los cuales se denominan cuadrados medios corregidos de la regresión y del error respectivamente (\(MS_{R}\) y \(MS_{Res}\))

A continuación una tabla que resume los hallazgos más importantes de todo esto, pero generalizada para el caso de regresión múltiple.

Además considerando que \(E(MS_{Res})={\sigma^2}_e\) y \(E(MS_{R})={\sigma^2}_e + {\beta_1}^2S_{xx}\)

Podríamos afirmar que si el valor \(F_0\) propuesto arriba llegase a ser muy grande, seguramente sería indicativo de que \(\beta_1 \ne 0\). Igualmente si \(\beta_1\) llegase a ser diferente de 0, entonces \(F_0\) se distribuye F de Snedecor con parámetro de no centralidad igual a:

\[\lambda = \frac{{\beta_1}^2S_{xx}}{{\sigma^2}_e}\]

Por consiguiente para probar la hipótesis \(H_0: \beta_1=0\) vs \(H_a: \beta_1 \ne0\), será suficiente con comparar \(F_0\) contra \(F_{\alpha,1,n-2}\), donde \(\alpha\) es el nivel de significancia para la prueba y por lo tanto se rechaza \(H_0\) si \(F_0>F_{\alpha,1,n-2}\)